The libraries
library(survival)
library(FRESA.CAD)
## Loading required package: Rcpp
## Loading required package: stringr
## Loading required package: miscTools
## Loading required package: Hmisc
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
##
## format.pval, units
## Loading required package: pROC
## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
library(Biobase)
## Loading required package: BiocGenerics
##
## Attaching package: 'BiocGenerics'
## The following object is masked from 'package:pROC':
##
## var
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, basename, cbind, colnames,
## dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
## grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
## rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
## union, unique, unsplit, which.max, which.min
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
##
## Attaching package: 'Biobase'
## The following object is masked from 'package:Hmisc':
##
## contents
## The following object is masked from 'package:miscTools':
##
## rowMedians
op <- par(no.readonly = TRUE)
pander::panderOptions('digits', 3)
pander::panderOptions('table.split.table', 400)
pander::panderOptions('keep.trailing.zeros',TRUE)
if (!require("BiocManager", quietly = TRUE))
{
install.packages("BiocManager")
BiocManager::install("seventyGeneData")
}
## Bioconductor version '3.15' is out-of-date; the current release version '3.17'
## is available with R version '4.3'; see https://bioconductor.org/install
library(seventyGeneData)
data(vanDeVijver)
class(vanDeVijver)
## [1] "ExpressionSet"
## attr(,"package")
## [1] "Biobase"
RRPlot Cox Model
timeinterval <- 5 # Five years
h0 <- sum(pdata$TTMevent & pdata$RFS <= timeinterval)
h0 <- h0/sum((pdata$RFS > timeinterval) | (pdata$TTMevent==1))
mcox <- coxph(Surv(RFS,TTMevent)~C1used,pdata)
pander::pander(summary(mcox)$coefficients)
| C1used |
-1.5 |
0.224 |
0.263 |
-5.69 |
1.3e-08 |
index <- predict(mcox,pdata)
rdata <- cbind(pdata$TTMevent,ppoisGzero(index,h0))
RRAnalysisCox <- RRPlot(rdata,atRate=c(0.90,0.95),
timetoEvent=pdata$RFS,
title="NIK: Breast Cancer",
ysurvlim=c(0.00,1.0),
riskTimeInterval=timeinterval)






par(op)
By Risk Categories
obsexp <- RRAnalysisCox$OERatio$atThrEstimates
pander::pander(obsexp)
| Total |
101 |
82.27 |
122.7 |
125.3 |
0.0285 |
| low |
6 |
2.20 |
13.1 |
14.7 |
0.0181 |
| 10% |
5 |
1.62 |
11.7 |
9.4 |
0.1893 |
| 5% |
90 |
72.37 |
110.6 |
99.0 |
0.3928 |
maxx <- 1.1*max(c(obsexp$Observed,obsexp$Expected))
minx <- 0.9*min(c(obsexp$Observed,obsexp$Expected))
plot(obsexp$Expected,obsexp$Observed,
xlim=c(minx,maxx),
ylim=c(minx,maxx),
main="Expected vs Observed",
ylab="Observed",
xlab="Expected",
col=rainbow(nrow(obsexp)),
log="xy")
errbar(obsexp$Expected,obsexp$Observed,obsexp$L.CI,obsexp$H.CI,add=TRUE,pch=0,errbar.col=rainbow(nrow(obsexp)),cex=0.75)
lines(x=c(1,maxx),y=c(1,maxx),lty=2)
text(obsexp$Expected,obsexp$Observed,rownames(obsexp),pos=2,cex=0.75)

RRPlot Cox Adjusted Model
This time we will include Lymph node status from pathology report and
Estrogen receptor alpha expression measurement from microarray
mcox <- coxph(Surv(RFS,TTMevent)~C1used*(ESR1 + Posnodes),pdata)
pander::pander(summary(mcox)$coefficients)
| C1used |
-0.403 |
0.668 |
0.629 |
-0.640 |
0.52186 |
| ESR1 |
0.123 |
1.131 |
0.255 |
0.481 |
0.63079 |
| Posnodesy |
-0.305 |
0.737 |
0.217 |
-1.401 |
0.16112 |
| C1used:ESR1 |
-1.913 |
0.148 |
0.739 |
-2.588 |
0.00966 |
| C1used:Posnodesy |
0.378 |
1.460 |
0.583 |
0.649 |
0.51661 |
index <- predict(mcox,pdata)
rdata <- cbind(pdata$TTMevent,ppoisGzero(index,h0))
RRAnalysisAdCox <- RRPlot(rdata,atRate=c(0.90,0.95),
timetoEvent=pdata$RFS,
title="Adjusted: Breast Cancer",
ysurvlim=c(0.00,1.0),
riskTimeInterval=timeinterval)






par(op)
By Risk Categories
obsexp <- RRAnalysisAdCox$OERatio$atThrEstimates
pander::pander(obsexp)
| Total |
101 |
82.27 |
122.7 |
122.43 |
0.0519 |
| low |
6 |
2.20 |
13.1 |
12.68 |
0.0658 |
| 10% |
5 |
1.62 |
11.7 |
7.37 |
0.5762 |
| 5% |
90 |
72.37 |
110.6 |
100.46 |
0.3183 |
maxx <- 1.1*max(c(obsexp$Observed,obsexp$Expected))
minx <- 0.9*min(c(obsexp$Observed,obsexp$Expected))
plot(obsexp$Expected,obsexp$Observed,
xlim=c(minx,maxx),
ylim=c(minx,maxx),
main="Ad. Expected vs Observed",
ylab="Observed",
xlab="Expected",
col=rainbow(nrow(obsexp)),
log="xy")
errbar(obsexp$Expected,obsexp$Observed,obsexp$L.CI,obsexp$H.CI,add=TRUE,pch=0,errbar.col=rainbow(nrow(obsexp)),cex=0.75)
lines(x=c(1,maxx),y=c(1,maxx),lty=2)
text(obsexp$Expected,obsexp$Observed,rownames(obsexp),pos=2,cex=0.75)

Calibrating the index
calprob <- CoxRiskCalibration(mcox,pdata,"TTMevent","RFS")
( 7.470999 , 626.7032 , 1.025293 , 101 , 118.8774 )
pander::pander(c(h0=calprob$h0,
Gain=calprob$hazardGain,
DeltaTime=calprob$timeInterval),
caption="Cox Calibration Parameters")
h0 <- calprob$h0
timeinterval <- calprob$timeInterval;
rdata <- cbind(pdata$TTMevent,calprob$prob)
RRAnalysisCalAdCox <- RRPlot(rdata,atRate=c(0.80,0.90),
timetoEvent=pdata$RFS,
title="Cal. NIK: Breast Cancer",
ysurvlim=c(0.00,1.0),
riskTimeInterval=timeinterval)






par(op)
By Risk Categories
obsexp <- RRAnalysisCalAdCox$OERatio$atThrEstimates
pander::pander(obsexp)
| Total |
101 |
82.27 |
122.7 |
120.1 |
0.0829 |
| low |
11 |
5.49 |
19.7 |
19.7 |
0.0539 |
| 20% |
10 |
4.80 |
18.4 |
13.6 |
0.4145 |
| 10% |
80 |
63.44 |
99.6 |
84.9 |
0.6639 |
maxx <- 1.1*max(c(obsexp$Observed,obsexp$Expected))
minx <- 0.9*min(c(obsexp$Observed,obsexp$Expected))
plot(obsexp$Expected,obsexp$Observed,
xlim=c(minx,maxx),
ylim=c(minx,maxx),
main="Cal. Expected vs Observed",
ylab="Observed",
xlab="Expected",
col=rainbow(nrow(obsexp)),
log="xy")
errbar(obsexp$Expected,obsexp$Observed,obsexp$L.CI,obsexp$H.CI,add=TRUE,pch=0,errbar.col=rainbow(nrow(obsexp)),cex=0.75)
lines(x=c(1,maxx),y=c(1,maxx),lty=2)
text(obsexp$Expected,obsexp$Observed,rownames(obsexp),pos=2,cex=0.75)

Comparing Risks
Comparing concordance Index
## Correlation Index
cindex <- RRAnalysisCI$c.index$cstatCI
## Cox Index
cindex <- rbind(cindex,RRAnalysisCox$c.index$cstatCI)
## Adjusted Cox Index
cindex <- rbind(cindex,RRAnalysisAdCox$c.index$cstatCI)
## Adjusted and Calibrated Cox Index
cindex <- rbind(cindex,RRAnalysisCalAdCox$c.index$cstatCI)
rownames(cindex) <- c("CI","Cox","Adj. Cox","Cal. Adj. Cox")
pander::pander(cindex)
| CI |
0.698 |
0.698 |
0.650 |
0.743 |
| Cox |
0.698 |
0.698 |
0.649 |
0.744 |
| Adj. Cox |
0.707 |
0.706 |
0.658 |
0.750 |
| Cal. Adj. Cox |
0.707 |
0.707 |
0.664 |
0.749 |
Comparing Risk Ratios Index
## Correlation Index
RRratio <- c(RR=RRAnalysisCI$keyPoints$RR[1],
LCI=RRAnalysisCI$keyPoints$RR_LCI[1],
UCI=RRAnalysisCI$keyPoints$RR_UCI[1])
## Cox Index
RRratio <- rbind(RRratio,c(RR=RRAnalysisCox$keyPoints$RR[1],
LCI=RRAnalysisCox$keyPoints$RR_LCI[1],
UCI=RRAnalysisCox$keyPoints$RR_UCI[1]))
## Adjusted Cox Index
RRratio <- rbind(RRratio,c(RR=RRAnalysisAdCox$keyPoints$RR[1],
LCI=RRAnalysisAdCox$keyPoints$RR_LCI[1],
UCI=RRAnalysisAdCox$keyPoints$RR_UCI[1]))
## Adjusted and Calibrated Cox Index
RRratio <- rbind(RRratio,c(RR=RRAnalysisCalAdCox$keyPoints$RR[1],
LCI=RRAnalysisCalAdCox$keyPoints$RR_LCI[1],
UCI=RRAnalysisCalAdCox$keyPoints$RR_UCI[1]))
rownames(RRratio) <- c("CI","Cox","Adj. Cox","Cal. Adj. Cox")
pander::pander(RRratio)
| CI |
3.81 |
2.08 |
6.96 |
| Cox |
3.81 |
2.08 |
6.96 |
| Adj. Cox |
3.42 |
1.93 |
6.07 |
| Cal. Adj. Cox |
3.24 |
2.10 |
4.98 |
Comparing logRank values
## Correlation Index
SurvDif <- c(chisq=RRAnalysisCI$surdif$chisq,pvalue=RRAnalysisCI$surdif$pvalue)
## Cox Index
SurvDif <- rbind(SurvDif,c(chisq=RRAnalysisCox$surdif$chisq,pvalue=RRAnalysisCox$surdif$pvalue))
## Adjusted Cox Index
SurvDif <- rbind(SurvDif,c(chisq=RRAnalysisAdCox$surdif$chisq,pvalue=RRAnalysisAdCox$surdif$pvalue))
## Adjusted and Calibrated Cox Index
SurvDif <- rbind(SurvDif,c(chisq=RRAnalysisCalAdCox$surdif$chisq,pvalue=RRAnalysisCalAdCox$surdif$pvalue))
rownames(SurvDif) <- c("CI","Cox","Adj. Cox","Cal. Adj. Cox")
pander::pander(SurvDif)
| CI |
28.1 |
7.97e-07 |
| Cox |
28.1 |
7.97e-07 |
| Adj. Cox |
28.5 |
6.46e-07 |
| Cal. Adj. Cox |
43.5 |
3.65e-10 |
Comparing Sensitivity
## Correlation Index
sensi <- RRAnalysisCI$ROCAnalysis$sensitivity
## Cox Index
sensi <- rbind(sensi,RRAnalysisCox$ROCAnalysis$sensitivity)
## Adjusted Cox Index
sensi <- rbind(sensi,RRAnalysisAdCox$ROCAnalysis$sensitivity)
## Adjusted and Calibrated Cox Index
sensi <- rbind(sensi,RRAnalysisCalAdCox$ROCAnalysis$sensitivity)
rownames(sensi) <- c("CI","Cox","Adj. Cox","Cal. Adj. Cox")
pander::pander(sensi)
| CI |
0.891 |
0.813 |
0.944 |
| Cox |
0.891 |
0.813 |
0.944 |
| Adj. Cox |
0.891 |
0.813 |
0.944 |
| Cal. Adj. Cox |
0.792 |
0.700 |
0.866 |
Comparing Specificity
## Correlation Index
speci <- RRAnalysisCI$ROCAnalysis$specificity
## Cox Index
speci <- rbind(speci,RRAnalysisCox$ROCAnalysis$specificity)
## Adjusted Cox Index
speci <- rbind(speci,RRAnalysisAdCox$ROCAnalysis$specificity)
## Adjusted and Calibrated Cox Index
speci <- rbind(speci,RRAnalysisCalAdCox$ROCAnalysis$specificity)
rownames(speci) <- c("CI","Cox","Adj. Cox","Cal. Adj. Cox")
pander::pander(speci)
| CI |
0.397 |
0.328 |
0.469 |
| Cox |
0.397 |
0.328 |
0.469 |
| Adj. Cox |
0.397 |
0.328 |
0.469 |
| Cal. Adj. Cox |
0.572 |
0.499 |
0.643 |
Comparing O/E
OERatio <- NULL
## Cox Index
OERatio <- rbind(OERatio,RRAnalysisCox$OERatio$estimate)
## Adjusted Cox Index
OERatio <- rbind(OERatio,RRAnalysisAdCox$OERatio$estimate)
## Adjusted and Calibrated Cox Index
OERatio <- rbind(OERatio,RRAnalysisCalAdCox$OERatio$estimate)
rownames(OERatio) <- c("Cox","Adj. Cox","Cal. Adj. Cox")
pander::pander(OERatio)
| Cox |
0.806 |
0.656 |
0.979 |
0.0285 |
| Adj. Cox |
0.825 |
0.672 |
1.002 |
0.0519 |
| Cal. Adj. Cox |
0.841 |
0.685 |
1.022 |
0.0829 |
Comparing O/Acum
OARatio <- NULL
## Cox Index
OARatio <- rbind(OARatio,RRAnalysisCox$OARatio$estimate)
## Adjusted Cox Index
OARatio <- rbind(OARatio,RRAnalysisAdCox$OARatio$estimate)
## Adjusted and Calibrated Cox Index
OARatio <- rbind(OARatio,RRAnalysisCalAdCox$OARatio$estimate)
rownames(OARatio) <- c("Cox","Adj. Cox","Cal. Adj. Cox")
pander::pander(OARatio)
| Cox |
1.15 |
0.936 |
1.40 |
0.165 |
| Adj. Cox |
1.18 |
0.961 |
1.43 |
0.104 |
| Cal. Adj. Cox |
1.01 |
0.823 |
1.23 |
0.881 |
Comparing NetBenefit
NetBen <- NULL
## Cox Index
NetBen <- rbind(NetBen,RRAnalysisCox$keyPoints$NetBenefit)
## Adjusted Cox Index
NetBen <- rbind(NetBen,RRAnalysisAdCox$keyPoints$NetBenefit)
## Adjusted and Calibrated Cox Index
NetBen <- rbind(NetBen,RRAnalysisCalAdCox$keyPoints$NetBenefit)
colnames(NetBen) <- rownames(RRAnalysisCox$keyPoints)
rownames(NetBen) <- c("Cox","Adj. Cox","Cal. Adj. Cox")
pander::pander(NetBen)
| Cox |
0.223 |
0.232 |
0.222 |
0.239 |
0.252 |
-0.02795 |
| Adj. Cox |
0.230 |
0.249 |
0.202 |
0.252 |
0.274 |
0.00617 |
| Cal. Adj. Cox |
0.162 |
0.190 |
0.162 |
0.215 |
0.239 |
0.01354 |